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Abstract 

The problem of time synchronization in dense wireless networks is considered. Well established synchro¬ 
nization techniques suffer from an inherent scalability problem in that synchronization errors grow with 
an increasing number of hops across the network. In this work, a model for communication in wireless 
networks is first developed, and then the model is used to define a new time synchronization mechanism. 
A salient feature of the proposed method is that, in the regime of asymptotically dense networks, it can 
average out all random errors and maintain global synchronization in the sense that all nodes in the 
multi-hop network can see identical timing signals. This is irrespective of the distance separating any 
two nodes. 


Index Terms 

Cooperation in networks, large network asymptotics, relay networks, scalability, sensor networks, time synchro¬ 
nization, wireless communications. 

I. Introduction 

A. Time Synchronization in Large Distributed Systems 

The problem of time synchronization in large distributed systems consists of giving all the physically 
disjoint elements of the system a common time scale on which to operate. This common time scale is 
usually achieved by periodically synchronizing the clock at each element to a reference time source, so 
that the local time seen by each element of the system is approximately the same. Time synchronization 
plays an important role in many systems in that it allows the entire system to cooperate and function as 
a cohesive group. 
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Time synchronization is an old problem [26], but the question of scalability is not. Recent advances in 
sensor networks show a clear trend towards the development of large scale networks with high node 
density. For example, a hardware simulation-and-deployment platform for wireless sensor networks 
capable of simulating networks with on the order of 100,000 nodes was recently developed [24]. As 
well, for many years the Smart Dust project sought to build cubic-millimeter motes for a wide range of 
applications [43]. Also, there is work in progress on the drastic miniaturization of power sources [27]. 
These developments (and many others) indicate that large scale, high density networks are on the horizon. 

Large scale, high density networks have applications in a variety of situations. Consider, for example, 
the military application of sniper localization. Large numbers of wireless nodes can be deployed to find 
the shooter location as well as the trajectory of the projectile [1]. Since the effective range of a long-range 
sniper rifle can be nearly 2km, in order to fully track the trajectory of the projectile it may be essential 
that our deployed network be tightly synchronized over distances of a few kilometers. Another example 
might be the implementation of a distributed radio for communication. In extracting information from a 
deployed sensor network, it may be beneficial for the nodes to cooperatively transmit information to a 
far away receiver [6], [7], [20]. Such an application would require that nodes across the network be well 
synchronized. As a result, a need for the synchronization of large distributed systems is very real and 
one that requires careful study to understand the fundamental performance limits on synchronization. 

B. Approaches to Synchronization and the Limitations 

The synchronization of large networks has been studied in fields ranging from biology to electrical 
engineering. The study of synchronous behavior has generally taken one of two approaches. The first 
approach is to consider synchronization as an emergent behavior in complex networks of oscillators. 
In that work, models are developed to describe natural phenomena and synchronization emerges from 
these models. The second approach is to develop and analyze algorithms that synchronize engineering 
networks. Nodes are programmed with algorithms that estimate clock skew and clock offset to achieve 
network synchronization. However, both of these approaches have significant limitations. 

1) The Emergence of Synchronous Behavior: Emergent synchronization properties in large populations 
has been the object of intense study in the applied mathematics ([30], [41]), physics ([3], [4], [5], [9], 
[12], [14], [16], [25]), and neural networks ([21], [37]) literature. These studies were motivated by a 
number of examples observed in nature: 

• In certain parts of south-east Asia, thousands of male fireflies congregate in trees and flash in 
synchrony at night [2]. 
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• Pacemaker cells of the heart, which on average cause 80 contractions a minute during a person’s 
lifetime [22]. 

• The insulin-secreting cells of the pancreas [35]. 

For further information and examples, see [32], [40], [31], [42], and the references therein. 


A number of models have been proposed to explain the emergence of synchrony, but perhaps one of the 
most successful and well known is the model of pulse-coupled oscillators by Mirollo and Strogatz [32], 


based on dynamical systems theory. Consider a function / : [0,1] —> [0,1] that is smooth, monotone 
increasing, concave down (i.e., /' > 0 and /" < 0), and is such that /(O) = 0 and /(I) = 1. Consider 
also a phase variable cj) such that d(j)/dt = where T is the period of a cycle. Then, each element in 
a group of N oscillators is described by a state variable Xi G [0,1] and a phase variable (pi G [0,1] as 
follows: 

• In isolation, Xi{t) = 

• If pi(t) = 0 then Xi{t) = 0, and if (pi(t) = 1 then Xi(t) = 1. 

• When 3:i(fo) = 1 for any of the Ts and some time ^o^ then for all other 



f ^{xj{pj(to)) + Si), Xj{pj{to)) +Si<l 
1, Xj{pj{to)) +Si> 1, 


where fg denotes an infinitesimal amount of time after to- That is, oscillator i reaching the end of 
a cycle causes the state of all other oscillators to increase by the amount Si, and the phase variable 
to change accordingly. 

The state variable Xi can be thought of as a voltage. Charge is accumulated over time according to 
the nonlinearity / and it discharges once it reaches full charge, resetting the charging process. Upon 
discharging, it causes all other charges to increase by a fixed amounf of Si, up fo fhe discharge poinf. For 
fhis model, if is proved in [32] fhaf for all N and for almosf all initial conditions, fhe system evenfually 
becomes synchronized. 

For fhe nefwork fo converge info a synchronous sfafe, one key assumption is fhaf fhe behavior of every 
single oscillator is governed by fhe same function /(•). This means fhaf all oscillators musf have fhe same 
frequency. From fhe liferafure, if appears fhaf fhis requiremenf is nearly always needed. As far as we are 
aware, for a fully synchronous behavior fo emerge, fhe oscillafors need to have fhe same, or nearly fhe 
same, oscillation frequencies. 

The need for nearly identical oscillators presenfs a significanl limifafion for emergen! synchronizafion. 
This emergence of synchrony is clearly desirable and if has been considered for communicafion and sensor 
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networks in [17], [18], [28]. However, whether or not these techniques can he adapted to synchronize 
networks with nodes that have arbitrary oscillator frequencies (clock skew) is still unclear. Thus, in 
order to overcome this limitation and find techniques capable of synchronizing a more general class of 
networks, we turn to algorithms designed to estimate certain unknown parameters such as clock skew. 

2) Estimation of Synchronization Parameters and the Scalability Problem: There have been many 
synchronization techniques proposed for use in sensor networks. These algorithms generally allow each 
node to estimate its clock skew and clock offset relative to the reference clock. Reference Broadcast 
Synchronization (RBS) [8] eliminates transmitter side uncertainties by having a transmitter broadcast 
reference packets to the surrounding nodes. The receiving nodes then synchronize to each other using 
the arrival of the reference packets as synchronization events. Tiny-Sync/Mini-Sync [36] and the Timing- 
sync Protocol for Sensor Networks (TPSN) [11] organize the network into a hierarchial structure and 
the nodes are synchronized using pair-wise synchronization. In lightweight tree-based synchronization 
(LTS) [13], pair-wise synchronization is also employed but the goal of LTS is to reduce communication 
and computation requirements by taking advantage of relaxed accuracy constraints. The Flooding Time 
Synchronization Protocol (FTSP) [29] achieves one-hop synchronization by having a root node broadcast 
timing information to surrounding nodes. These surrounding nodes then proceed to broadcast their 
synchronized timing information to nodes beyond the broadcast domain of the root node. This process 
can continue for multi-hop networks. 

The problem with each of these traditional synchronization techniques is that synchronization error 
will increase with each hop. Since each node is estimating certain synchronization parameters, i.e. clock 
skew, there will be inherent errors in the estimate. As a result, a node multiple hops away from the node 
with the reference clock will be estimating its parameters from intermediate nodes that already have 
estimation errors. Therefore, this introduces a fundamental scalability problem: as the number of hops 
across the network grows, the synchronization error across the network will also grow. 

Current trends in network technology are clearly moving us in the direction of large, multi-hop 
networks. First, sensors are decreasing in size and this size decrease will most likely be accompanied by a 
decrease in communication range. Thus, more hops will be required to traverse a network deployed over 
a given area. Second, as we deploy networks over larger and larger areas, then for a given communication 
range, the number hops across the network will also increase. In either case, the increased number of 
hops required to communicate across the network will increase synchronization error. Therefore, it is 
essential that we develop techniques than can mitigate the accumulation of synchronization error over 
multiple hops. 
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C. Spatial Averaging and Synchronization 

1) Cooperation through Spatial Averaging: To decrease the error increase in each hop, we need to 
decrease the estimation error. There are two primary ways of achieving this. First, each node can increase 
the amount of timing information it obtains from neighboring nodes. For example, from a received timing 
packet, the node may be able to construct a data point telling it the approximate time at the reference 
clock and the corresponding time at its local clock. Using a collection of these data points, the node can 
estimate clock skew and clock offset. So instead of using, say, five packets with timing information, a 
node can wait for ten packets. More data points will generally give better estimates. The drawback to 
such an approach is the increase in the number of packet exchanges. 

The second way in which to reduce estimation error is to increase the quality of each data point 
obtained by the nodes. This can be achieved through improving packet exchange algorithms and time 
stamping techniques. However, we believe that there is one fundamentally new approach to improving 
data point quality that has not be carefully studied. This is to use spatial averaging to improve the quality 
of each data point. 

The motivation for this approach is very simple. Assume that each node has many neighbors. If all 
nodes in the network are to be synchronized, then the neighbors of any given node will also have 
synchronization information. Is it possible to simultaneously use information from all the neighbors to 
improve the quality of a timing observation made by a node? Furthermore, it would seem to make 
sense that with more neighbors, hence more available timing information, the quality of the constructed 
data point should improve. If this is indeed the case, then achieving synchronization through the use of 
spatial averaging will provide a fundamentally new trade-off in improving synchronization performance. 
Network designers would simply be able to increase the number and density of nodes to obtain better 
network synchronization. The study of cooperative time synchronization using spatial averaging is the 
focus of this work. 

2) Model and Technique: To obtain a model for developing cooperative synchronization in large 
wireless networks, we begin by looking at the signals observed by a node in a network with N nodes 
uniformly deployed over a fixed finite area. To start, we assume propagation delay to negligible (the 
general case is considered in Section V). All nodes transmit a pulse p{t) and a node j will see a signal 
Aj^N{t) which is the superposition of all these pulses, 

N . „ 

i=l 
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In this expression, p{t) is the basic pulse transmitted hy each node (assumed to he the same for all 
nodes), tq is the ideal pulse transmit time, hut since we assume imperfect time synchronization among 
the nodes we have Tj modelling random errors in the pulse transmission time. Kj^i models the amplitude 
loss in the signal transmitted hy the zth node. Amax is the maximum magnitude transmitted hy a node. 
We scale each node’s transmission hy N so that as the network density grows, the total power radiated 
does not grow unbounded. This model thus describes the received signal seen at a node j for a network 
with N nodes and this holds for any N. Increasing N will have two effects: (a) node density will increase 
since the network area is fixed and (b) node signal transmission magnitude will decrease due to the 1/N 
scaling. Therefore, by increasing N this model allows us to study the scalability of networks as node 
density grows and node size decreases. 

Given that these are the signals observed at each node, we ask: is it possible for Aj^]\f(t) to encode 
a time synchronization signal that will enable all nodes in the network to synchronize their clocks with 
bounded error, as N ^ oo? The answer is yes, and the key to proving all our results is the law of large 
numbers. 

Our key idea is the following. If all nodes were able to determine when time tq (in the reference 
time) arrives, then by transmitting p{t) at time tq, the signal observed at any node j would be p{t — 
''"o) , which is a suitably scaled version of p{t) centered at tq. In reality however, there 

will be some error in the determination of tq, which we account for by allowing for a node-dependent 
random error Tj. But, if the distribution of Tj satisfies certain conditions, then the effects of that timing 
error can be averaged out. A pictorial representation of why this should be the case is shown in Fig. 1. 



Fig. 1. Assume N square waves are transmitted (one by each node) at random times. These times have the properties 
that they all have the same mean, a small variance compared to the duration of the wave, and their distribution is 
symmetric. Then, under the assumption of N large, it follows from the Law of Large Numbers that the observed 
signal is going to be a smoothed version of the square wave, in which the center zero-crossing will correspond to 
the location of the mean of the random times. 
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Therefore, intuitively we can see how the technique of cooperative time synchronization using spatial 
averaging can average out the inherent timing errors in each node. Even though there is randomness 
and uncertainty in each node’s estimates, hy using cooperation among a large number of nodes it is 
possible to recover deterministic parameters from the resulting aggregate waveform (such as the location 
of certain zero-crossings) in the limit as node density grows unbounded. Thus more nodes will give us 
better estimates. This is because the random waveform converges to a deterministic one as more and 
more nodes cooperatively generate an aggregate waveform. At the same time, the average power required 
by each node will decrease since smaller nodes send smaller signals. Therefore, by programming suitable 
dynamics into the nodes, in this paper we show how it is possible to generate an aggregate output signal 
with equispaced zero-crossings in the limit of asymptotically dense networks. Thus, the detection of 
these zero-crossings plays the same role as that of an externally generated time reference signal based 
on which all nodes can synchronize. 

We develop this synchronization technique in three main steps. One, we set up the model for Aj w(0- 
Two, we specify characteristics of the model (i.e. the distribution of Tj) that allow us to prove desirable 
properties of the aggregate waveform (such as a center zero-crossing at tq). Three, we develop the 
estimators needed for our synchronization technique and show that the estimators give us the desired 
characteristics. 

D. Main Contributions and Organization of the Paper 

The main contributions presented in this paper are the following; 

• The definition of a probabilistic model for the study of the time synchronization problem in wireless 
networks. This model does contain the classical Mirollo-Strogatz model as a special case, but its 
formulation and the tools used to prove convergence results are of a completely different nature 
(purely probabilistic, instead of based on the theory of dynamical systems). 

• Using this model, we provide a rigorous analysis of a new cooperative time synchronization technique 
that employs spatial averaging and has favorable scaling properties. As the density of nodes increases, 
synchronization performance improves. In particular, in the limit of infinite density, deterministic 
parameters for synchronization can be recovered. 

• We show that cooperative time synchronization works perfectly for negligible propagation delay. 
When propagation delay is considered, we find thaf asymmetries at the boundaries reveal some 
limitations that need to be carefully considered in designing algorithms that take advantage of spatial 
averaging. 
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In analyzing the proposed cooperative time synchronization technique, our goal is to show that the pro¬ 
posed technique can average out all random error and provide deterministic parameters for synchronization 
as node density grows unbounded. This asymptotic result can he viewed as a convergence in scale to 
synchrony. The result serves as a theoretical foundation for allowing a new trade-off between node density 
and synchronization performance. In particular, higher node density can yield better synchronization. 

The rest of this paper is organized as follows. The general model is presented in Section II. Of particular 
interest here is Section II-E, where we show how our model contains the model of Mirollo and Strogatz 
for pulse-coupled oscillators as a special case [32]. In Section III we specialize the general model for 
our synchronization setup and develop waveform properties that will be used in time synchronization. 
In Section IV we develop the cooperative time synchronization technique for no propagation delay. We 
extend the cooperative synchronization ideas to the case of propagation delay in Section V. The paper 
concludes in Section VI with a detailed discussion on the scalability issue and how the technique proposed 
in this work lays the theoretical foundation for a general class of cooperative time synchronization 
techniques that use spatial averaging. 


II. System Model 

A. Clock Model 

We consider a network with N nodes uniformly distributed over a fixed finite area. The behavior of 
each node i is governed by a clock Cj that counts up from 0. The introduction of q is important since 
it provides a consistent timescale for node i. By maintaining a table of pulse-arrival times, node i can 
utilize the arrival times of many pulses over an extended period of time. 

The clock of one particular node in the network will serve as the reference time and to this clock we 
wish to synchronize all other nodes. We will call the node with the reference clock node 1 and the clocks 
of other nodes are defined relative to the clock of node 1. We never adjust the frequency or offset of the 
local clock Ci because we wish to maintain a consistent time scale for node i. 

The clock of node 1, ci, will be defined as ci{t) = t where t G [0, oo). Taking ci to be the reference 
clock, we now define the clock of any other arbitrary node i, Ci. We define Cj as 

Ci{t) = ai{t - Ki)+ ( 1 ) 

where 

• Aj is an unknown offset between the start times of q and ci. 
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• ttj > 0 is a constant and for each i, at G [aiow, Oup] where aup, aiow > 0 are finite. This hound on 
Ui means that the frequency offsets between any two nodes can not he arbitrarily large. 

• 'l'i(t) is a stochastic process modeling random timing jitter. 

Thus, this model assumes that there is a bounded constant frequency offset between the oscillators of 
any two nodes as well as some random clock jitter. 

It is important to note that node 1 does not have to be special in any way; its clock is simply a reference 
time on which to define the clocks of the other nodes. This means that our clock model actually describes 
the relative relationship of all the clocks in the network by using an arbitrary node’s clock as a reference. 

B. Pathloss Only Model 

1) A Random Model for Pathloss: From Section I-C.2, we see that we are interested in studying the 
aggregate waveform observed at a node j. As a result, we are only concerned with the aggregate signal 
magnitude and do not care about the particular signal contribution from each surrounding node. With 
this in mind, we can develop a random model for pathloss that, for dense networks, gives the appropriate 
aggregate signal magnitude at node j. Such a model is ideal for our situation since we are studying 
asymptotically dense networks. 

We start with a general pathloss model K{d), where Q < K{d) < 1 for all distances d > 0, is non¬ 
increasing and continuous. K{d) is a fraction of the transmitted magnitude seen at distance d from the 
transmitter. For example, if the receiver node j is at distance d from node i, and node i transmits a 
signal of magnitude A, then node j will hear a signal of magnitude AK{d). We derive K{d) from a 
power pathloss model since any pathloss model captures the average received power at a given distance 
from the transmitter. This average received power is perfect for modelling received signal magnitudes 
in our problem setup since we are considering asymptotically dense networks. Due to the large number 
of nodes at any given distance d from the receiver, using the average received magnitude at distance d 
as the contribution from each node at that distance will give a good modelling of the amplitude of the 
aggregate waveform. 

The random pathloss variable Kj will be derived from K{d). To understand how Kj and K{d) are 
related, we give an intuitive explanation of the meaning of Kj as follows: the Pr[Kj G {k,k + A)] 
is the fraction of nodes at distances d from node j such that K{d) G (fc, A: + A), where A is a small 
constant. This means that, roughly speaking, for any given scaling factor Kj = k, (A:) A is the fraction 
of received signals with magnitude scaled by approximately k, where fx^ {k) is the probability density 
function of Kj. Thus, if we scale the transmit magnitude A from every node i by an independent Kj, 
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then as the number of nodes, N, gets large, node j will see NfKj{k)A signals of approximate magnitude 
Ak, and this holds for all k in the range of Kj. This is because taking a large number of independent 
samples from a distribution results in a good approximation of the distribution. 

Thus, this intuition tells us that scaling the magnitude of the signal transmitted from every node i 
by an independent sample of the random variable Kj gives an aggregate signal at node j that is the 
same magnitude as if we generated the signal using K{d) directly. Even though the signals from two 
nodes at the same distance from a receiver have correlated magnitudes, we do not care about the signal 
magnitude from any particular node but only that an appropriate number of all possible received signal 
magnitudes contribute to the aggregate waveform. For a receiving node j, we choose therefore to work 
with the random variable Kj instead of directly with K{d) because, for the goals of this paper, doing so 
has two major advantages: (a) we can obtain desirable limit results by placing very minimal restrictions 
on the distribution of the Kj’s (and hence on K{d)) and (b) we can apply tools from probability theory 
(basically, the strong law of large numbers) to carry out our analysis. 

2) Definition of Kj: From the above intuition we can define the cumulative distribution function of 
Kj as 


FK,{k) = Fr{Kj<k) 


0 

AT-A{j,f) 

1 


k e (— oo, 0) 

1-^ te[o,i] 

k G (1, oo) 


( 2 ) 


where 

• At is the total area of the network, 

• A{j, a) is the area of the network contained in a circle of radius a centered at node j, 

• f = supjd : K{d) > k}. 

From the above discussion we see that the distribution of Kj is only a function of node j, the receiving 
node. We illustrate the relationship among node j, K(d), r, and Fk^ (k) in Fig 2. We sometimes write Kjj 
with i used to index each node surrounding node j. i is thus indexing a sequence of independent random 
variables Kjj for fixed j. Therefore, for a given j, Kj/s are independent and identically distributed 
(i.i.d.) with a cumulative distribution function given by (2) for all i. 

We assume that Kj has the following properties: 

• Kj is independent from 'k/(f) for all j, I, and t. 


• 0 < Kj < 1, 0 < E{Kj) < 1, and 'V3x{Kj) < 1. 


The requirements on the random variable Kj places restrictions on the model K{d). Any function K{d) 
that yields a Kj with the above requirements can be used to model pathloss. 
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Network Area , Pathloss Function 

d 



Cumulative Distribution Function of K j 

Fig. 2. An illustration of the cumulative distribution function Fk^ [k) is shown in the bottom-right figure. For a 
given scaling value k, Fk^ {k) is defined to be 1 — {A{j, r) jAx), where the relationship between f and k is shown 
in the top-right figure. The area A{j, f) and its relation to node j is shown in the top-left figure. 


C. Delay and Pathloss Model 

In this section we develop a more complex model to simultaneously model propagation delay and 
pathloss. This leads to the joint development of the delay random variable Dj and a corresponding 
pathloss random variable Kj. 

1) Correlation Between Delay and Pathloss: Since we want to develop a model for both pathloss 
and time delay, we start by keeping the pathloss function K{d) defined in Section II-B. The general 
delay model assumes a function 5{d) that models the time delay as a function of distance. 5[d) describes 
the time in terms of ci that it takes for a signal to propagate a distance d. For example, if node i and 
node j are distance do apart, then a pulse sent by node i at time ci = 0 will be seen at node j at 
time Cl = 6{do). We make the reasonable assumption that d{d) is continuous and strictly monotonically 
increasing for d > 0. 

As with the pathloss only model, we want to define a delay random variable Dj for each receiving 
node j. Recall fhaf this means that for every node j there is a random variable Dj associated with it 
since, in general, each node j will see different delays. There is a correlation between the delay random 
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variable Dj and the pathless random variable Kj. This correlation arises for two main reasons. First, 
since in Section II-B we define K{d) to be monotonically decreasing and continuous, it is possible for 
K{d) = 0 for d G [i?, oo), R > 0. This might be the case for a multi-hop network. In this situation, 
there will be a set of nodes whose transmissions will never reach node j (i.e. infinite delay) even though 
according to 6{d) these nodes should contribute a signal with finite delay. Second, a small Kj value 
would represent a signal from a far away node. As a result, the corresponding Dj value should be large 
to reflect large delay. Therefore, keeping these two points in mind, we proceed to develop a model for 
both pathloss and propagation delay. 

2) Definition of Dj and Kj: We define the cumulative distribution function of Dj as 


Fdj (x) = Pr(T>j < x) 


0 

A'j' 

a{x - 6iR)) + ^ 
1 


X G (—oo, 0) 

X G [0,S{R)] 

X G {6{R),6{R + AR)] 
X G {6{R + AR),oo) 


( 3 ) 


where r' = sup{r : d{r) < x}, AR > 0 is a constant, R = supjd : K(d) > 0}, and 

_ A{j,R) 

_ ^ At 

““ 6{R + AR) - 5{R) ■ 

Recall that A(j, a), defined in Section II-B, is the area of the network contained in a circle of radius a 
centered at node j and At is the total area of the network. Note that R can be infinite. 

Using the delay random variable Dj with the cumulative distribution function in (3), we define Kj as 


K,=K{6-HDj)), 


( 4 ) 


where K{-) is the deterministic pathloss function from Section II-B and 5“^ : [0,oo) —> [ 0 , 00 ) is the 
inverse function of the deterministic delay function (5(-). Note that exists since (5(-) is continuous 

and strictly monotonically increasing on [0,oo). 

3) Intuition Behind Dj and Kj: To understand the distribution of Dj, we need to consider the definition 
of Kj as well Recall that a signal arriving with delay Dj is scaled by the pathloss random variable Kj. 
Let us consider the cumulative distribution in two pieces, x G [0, and x G (<5(1?), 00 ). The case 

for X G (— 00 , 0 ) is trivial First, for x G [0, the probability that Dj takes a value less than or 

equal to x is simply the fraction of the network area around node j such that the nodes are at distances 
d with 5{d) < X. The intuition is the same as that for the development of Kj in Section II-B. Second, 
for X G {5{R),oo), the situation is more complex. Note that a transmitted signal from a node at distance 
d G {R,oo) from j will arrive at node j with infinite delay since K{d) =0 for d G {R,oo). Since any 
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delay values in x G {S{R), oo) correspond to distances d = G {R, oo), the corresponding scaling 

value will be zero because Kj and Dj are related by (4). As a result, it does not matter what delay values 
we assign to the fraction of the network area outside a circle of radius R centered at node j as long as 
their delay value x is such that G (ii, oo). Thus, we can arbitrarily choose a constant Aii value 

and construct a piecewise linear portion of the cumulative distribution function of Dj on x G {6{R),oo). 
The probability that Dj G {d{R), oo) will be the fraction of the network area outside a circle of radius R 
around node j. And since Dj G {d{R),oo) will have a corresponding Kj value that is zero, this fraction 
of nodes will not contribute to the aggregate waveform at node j. It is clear that the correlated Dj and Kj 
random variables work together to accurately model a signal arriving with both pathloss and propagation 
delay. An illustration of how K{d), 6{d), node j, and Fd^(x) are related can be found in Fig. 3. 



Fig. 3. From the top-left and bottom-left figures, we can see how K{d) determines the set of nodes surrounding 
node j that will contribute to the aggregate waveform at node j. This contributing set of nodes is related to Fjj. (x) 
through S{d) and this is illustrated in the top-right and bottom-right figures. 

We require that Dj is bounded, has finite expectation, and has finite variance for all j. Note that 
D,>0 by the requirement that d{d) > 0. As well, since the cumulative distribution in (3) is continuous, 
and often absolutely continuous, we assume that Dj has a probability density function /^^(x). When we 
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write Dj^i, the i indexes each node surrounding node j. Thus, the Dj/s are independent and identically 
distributed in i for a given j and have a cumulative distrihution given hy (3). Using the Kj and Dj 
developed in this section to simultaneously model pathloss and propagation delay, respectively, we will 
he able to closely approximate the received aggregate waveform at any node j as N ^ oo. 

To summarize, we see that our choice of the pathloss and delay random variables will depend on what 
we want to model. If we only consider pathloss and not propagation delay, then we will use the random 
variable Kj defined in Section II-B. If we account for both pathloss and delay, then we will use the delay 
random variable Dj in this section (Section II-C) and the pathloss random variable Kj defined by (4). 

D. Synchronization Pulses and the Pulse-Connection Function 

The exchange of pulses is the method through which the network will maintain time synchronization. 
Each node i will periodically transmit a scaled pulse Aip{t), where Ai is a constant and p{t), in general, 
can be any pulse. We call the interval of time during which a synchronization pulse is transmitted a 
synchronization phase. 

What each node does with a set of pulse arrival observations is determined by the pulse-connection 
function for node i. The pulse-connection function is a function that determines the time, in the 
time scale of q, when node i will send its nth pulse. It can be a function of the current value of Ci{t) 
and past pulse arrival times. This function basically determines how any node i reacts to the arrival of 
a pulse. 

E. An Example: Pulse-Coupled Oscillators 

The system model that we presented thus far is powerful because it is very general. In this section 
we show that it is a generalization of the pulse-coupled oscillator model proposed by Mirollo and 
Strogatz [32]. As a result, the results presented in that paper will hold under the simplified version 
of our model. 

1) Model Parameters for Pulse-Coupled Oscillators: In selling up the system model, Mirollo and 
Strogatz make four key assumptions: 

• Pathloss Model: The first assumption that is made is that there is all-to-all coupling among all N 
oscillators. This means that each oscillator’s transmission can be heard by all other oscillators. Thus, 
for our model we ignore pathloss, i.e. K{d) = 1, to allow any node’s transmission to be heard by 
each of the other N — 1 nodes. 
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• Delay Model: The second assumption is that there is instantaneous coupling. This assumption is the 
same as setting 6{d) = 0. In such a situation we would use our pathloss only model. 

• Synchronization Pulses: The third key assumption made in [32] is that there is non-uniform coupling, 
meaning that each of the N oscillators fire with strengths ei,... We modify the parameters in 
our model hy making node i transmit with magnitude Ai = e^. They also assume that any two pulses 
transmitted at different times will he seen hy an oscillator as two separate pulses. In our model, we 
may choose any pulse p{t) that has an arhitrarily short duration and each node will detect the pulse 
arrival time and pulse magnitude. 

• Clock Model: The fourth important assumption made hy Mirollo and Strogatz is that the oscillators 
are identical hut they start in arbitrary initial conditions. We simplify our clock model in (1) hy 
eliminating any timing jitter, i.e. 'l'i(f) = 0, and making the clocks identical hy setting = 1 for 
i = I,... ,N. We leave Aj in the model to account for the arbitrary initial conditions. We also 
assume that the phase variable in the pulse-coupled oscillator model increases at the same rate as 
our clock. That is, the time it takes the phase variable to go from zero to one and the time it takes 
our clock to count from one integer value to the next are the same. 

Now that we have identical system models, what remains is to modify our model to mimic the coupling 
action detailed in [32]. This is accomplished by defining a proper pulse-connection function 

2) Choice of Pulse-Connection Function: To match the coupling action in [32], we choose a pulse 
transmit time function Xf ..., fj that is a function of pulse receive times and 

also the time of node i’s (n — l)th pulse transmission time, z'^f^ is the time in terms of c, that node 
i receives its feth pulse since its last pulse transmission at In this case, X^^ will be a function 

that updates node Ts nth pulse transmission time each time node i receives a pulse. Let X^\{k) = 
^k-i ii ■ ■ ■ ’ ^Ti^^n -1 i) where it is node i’s nth pulse transmission time after observing k pulses 
since its last pulse transmission. Node i will transmit its pulse as soon as X^\ < Ci{t) where Ci{t) is 
node Ts current time. As soon as the node transmits a pulse at X^^ the function will reset and become 
X'fj^nff) = xW + 1. The node is now ready to receive pulses and at its first received pulse, the next 
transmission time will become i(^)- ^ni will thus be defined as 




-!)-[/ 


-1 


(ej + f{4\i - - 4k,i - ^ > 0 


+ 1 


( 5 ) 

( 6 ) 


where the pulse received at is a pulse of magnitude ej and the function / : [0,1] —> [0,1] is the 
smooth, monotonic increasing, and concave down function defined in [32]. 
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Equations (5) and (6) fundamentally say that each time node i receives a pulse, node i’s next trans¬ 
mission time will he adjusted. This is in line with the behavior of the coupling model descrihed hy 
Mirollo and Strogatz since each time an oscillator receives a pulse, its state variable is pulled up by e 
thus adjusting the time at which the oscillator will next fire. To see how equations (5) and (6) relate 
to the coupling model in [32], let us consider an example with two pulse coupled oscillators. Consider 
two oscillators A and B illustrated in Fig. 4. In Fig. 4(a), we have that oscillator A is at phase q and 



+1 


(a) 



x";(i) x°;(o) 

(b) 


Fig. 4. We illustrate the connection between the pulse-coupled oscillator coupling model and our clock model. In 
(a), oscillator B is just about to fire and oscillator A has phase q. In (b), oscillator B fires and increases the phase 
of oscillator A by d. This d increase in phase effectively decreases the time at which A will next fire. We capture 
this time decrease by decreasing the firing time of our node by an amount d. Thus, oscillator A and our node will 
fire at the same time. 


oscillator B is just about to fire. Below the pulse-coupled oscillator model we have a time axis for node i 
corresponding to our clock model going from time ^ to x'^_^ i + 1- Our time axis for node i models 
the behavior of oscillator A, that is, we want node i to behave in the same way as oscillator A under the 
influence of oscillator B. If oscillator B did not exist, then the phase variable q will match our clock in 
that q reaches 1 at the same time our clock reaches = x'^_^ ^ + 1 and oscillator A will fire af 

fhe same time our model fires. In Fig. 4(b), oscillator B has fired and has pulled the state variable of 
oscillator A up by e. This coupling has effectively pushed the phase of oscillator A to q + d and decreased 
the time before A fires. In fact, the time until oscillator A fires again is decreased by d. We can capfure 
fhis coupling in our model since we can calculate the lost time d. The time at which oscillator B fires is 
and it is clear that d = f~^{e + f{z\\ — x^_^ J) — (z'l^ — x^_^ ■). Thus, if the time that oscillator 
A will fire again is decreased by time d due to the pulse of B, then we adjust our node firing time by 
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decreasing the firing time to ^ + 1 — d. This is exactly the expression in (5) for A; = 1. 

This relationship between our model for calculating the node firing time and the pulse-coupled oscillator 
coupling model can he easily extended to N oscillators. 

We can see then that the pulse-coupled oscillator model proposed hy Mirollo and Strogatz in [32] is 
a special case of our model. Our model generalizes this pulse-coupled oscillator model hy considering 
timing jitter, pulses of finite width, propagation delay, non-identical clocks, and an ability to accommodate 
arbitrary coupling functions. 

III. Cooperative Time Synchronization Setup 

Just as we could specialize our model to the pulse-coupled oscillator model of Mirollo and Strogatz, 
we now specialize the model for our proposed synchronization technique. We start under the assumption 
of no propagation delay and develop the synchronization technique for this case. Propagation delay is 
considered in Section V. We proceed in three steps. In Section III-B, we specify the model for 
the received waveform at any node j. Second, in Section III-C, we prove that given certain characteristics 
of the model, has very useful limiting properties. Third, we show in Section IV that estimators 

(i.e., the pulse connection function) developed for our synchronization technique give the desired 

properties. 

A. System Parameters 

For our synchronization technique, we specialize the general model by making the following assump¬ 
tions on Ui and 'l'j(f) for f = 1... 

• A characterization of the {ui} is given by a known function fa{s) with s G [aiow^otup] that gives 
the percentage of nodes with any given a value. Thus, the fraction of nodes with a values in the 
range sq to si can be found by integrating /q(s) from sq to si. We assume that |/a(s)| < Ga, 
for some constant Ga- We keep this function constant as we increase the number of nodes in the 
network (N —> oo). Given any circle of radius R that intersects the network, the nodes within that 
circle will have a^’s that are characterized by /a(s). R is the maximum d such that K{d) > 0. This 
means that the set of nodes that any node j will hear from will have its ctj’s characterized by a 
known function. Note that R can be infinite, and in that case, any node j hears from all nodes in 
the network. Fundamentally, fa{s) means that as we increase node density, the new nodes have a 
parameters that are well distributed in a predictable manner. 
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• 'l'i(t) is a zero mean Gaussian process with samples 'l'i(to) ~ AA(0, cj^), for any to, and independent 
and identically distributed samples for any set of times [to,... ,1^], k a positive integer. We assume 
< oo and note that is defined in terms of the clock of node i. We assume that is Gaussian 
since the RMS (root mean square) jitter is characterized by the Gaussian distribution [34]. 

We maintain the full generality of the pathloss model from Section II-B. Note that throughout this work 
we assume no transmission delay or time-stamping error. This means that a pulse is transmitted at exactly 
the time the node intends to transmit it. We make this assumption since there will be no delay in message 
construction or access time [8] because our nodes broadcast the same simple pulse without worrying about 
collisions. Also, when a node receives a pulse it can determine its clock reading without delay since any 
time stamping error is small and can be absorbed into the random jitter. 

B. Signal Reception Model 

For our proposed synchronization technique, the aggregate waveform seen by node j at any time t is 

N . 

TN(«) = E - T.), (7) 

i=\ 

where A^>^{t) is the waveform seen at node j written in the time scale of ci and Ai = Amax/N for 
all i. Amax is the maximum transmit magnitude of a node. Tj is the random timing offset suffered by 
the zth node, which encompasses the random clock jitter and estimation error. This model says that 
each node i’s pulse transmission occurs at the ideal transmit time tq plus some random error Tt. In the 
next section. Section III-C, we find properties for Ti that will give us desirable properties in A'i^^{t). 
Then, in Section IV, we show that our proposed steady-state synchronization technique and its associated 
pulse-connection function will give us the desired properties. 

There are two comments about (7) that we want to make. First, note that even though we sum the 
transmissions from all N nodes in (7), we do not assume that node j can hear all nodes in the network. 
Recall from the pathloss model that if we have a multi-hop network, then there will be a nonzero 
probability that Kj^i = 0. Thus, node j will not hear from the nodes whose transmissions have zero 
magnitude. Second, it may be possible that the nodes are told that there are TV = vN nodes in the 
network while the actual number of functioning nodes is N. In which case, each node will transmit with 
signal magnitude Ai = Amax/{vN) and (7) will have a factor of 1/u. Other than for this factor, however, 
the theoretical results that follow are not affected. 

To model the quality of the reception of A^>^{t) by node j, we model the reception of a signal by 
defining a threshold 7. 7 is the received signal threshold required for nodes to perfectly resolve the pulse 
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arrival time. If the maximum received signal magnitude is less than 7 then the node does not make any 
observations and ignores the received signal waveform. We assume that 7 <C Amax- 


In our work we will assume that p{t) takes on the shape 



p{t) - < 0 t = 0,t< -Tnz, t > Tnz 

-q{-t) 0 <t <Tnz 


( 8 ) 


where Tnz > 0 is expressed in terms of ci. We assume q{t) > 0 for f e (—Tn^,0), q{t) 7 ^ 0 only on 
t G (—Tn 2 , 0 ), sup^|( 7 (f)| = 1, and q{t) is uniformly continuous on (—r„ 2 , 0 ). Thus, we see that p{t) 
has at most three jump discontinuities (at t = 0, —Tnz, Tnz)- Tnz should he chosen large compared to 
maxjfj?, i.e. af « Tnz, where af is the value of translated from the time scale of c* to ci. This 
way, over each synchronization phase, with high prohahility a zero-crossing will occur. For each node, 
the duration in terms of ci of a synchronization phase will he 2Tnz- Note that we assume Tm is a value 
that is constant in any consistent time scale. This means that even though nodes have different clocks, 
identical pulses are transmitted hy all nodes. We define a pulse to he transmitted at time t if the pulse 
makes a zero-crossing at time t. Similarly, we define fhe pulse receive (arrival) time for a node as fhe 
fime when fhe observed waveform firsf makes a zero-crossing. A zero-crossing is defined for signals fhaf 
have a posifive amplifude and fhen fransifion fo a negafive amplifude. If is fhe fime fhaf fhe signal firsf 


reaches zero. 


For fhe exchange of synchronizafion pulses, we assume fhaf nodes can fransmif pulses and receive 


signals at the same time. This simplifying assumption is not required for the ideas presented here to 


hold, but simplifies fhe presenfafion. We mention a way to relax this assumption in Section IV-D.l. 

In (7) and in the discussions above, we have focused on characterizing the aggregate waveform for 
any one synchronization phase. That is, (7) is the waveform seen by any node j for the synchronization 
phase centered around node I’s transmission at f = tq, tq a positive integer. We can, however, describe 
a synchronization pulse train in the following form. 



(9) 


N 


U=1 2=1 


where t^ is the integer value of t at the nth synchronization phase, and Tj „ is the error suffered by the ith 
node in the uth synchronization phase. We seek to create this pulse train with equispaced zero-crossings 
and use each zero-crossing as a synchronization event. An illustration of such a pulse train is shown in 
Fig. 5. For simplicity, however, most of the theoretical work is carried out on one synchronization phase. 
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time 


Fig. 5. An illustration of a pulse train with equispaced zero-crossings. The pulse at each integer value of t is 
an instance of Aj^oo{t) = limjv^oo V(^) three instances of Aj^ao{t) in the above figure with zero- 

crossings at f = 1,2, 3. We can control the zero-crossings of Aj oo(f) and choose to place it on an integer value of 
t. As a result, we can use these zero-crossings as synchronization events since they can be detected simultaneously 
by all nodes in the network. 





C. Desired Structural Properties of the Received Signal 

In this section, we characterize the properties of Tj that give us desirable properties in the aggregate 
waveform. From (7), the aggregate waveform seen at each node j in the network has the form 

1 ^ 

Awit) = ^ X] AmaxKip{t - To -Ti) (10) 

i=l 

We have dropped the j and ci for notational simplicity since in this section we deal solely with the 
received waveform at a node j in the time scale of ci- As we let the number of nodes grow unbounded 
(N —> oo), the properties of this limit waveform can be characterized by Theorem 1. These properties 
will be essential for asymptotic cooperative time synchronization. As a note, in Theorem 1 we present the 
case for Gaussian distributed Tj but similar results hold for arbitrary zero-mean, symmetrically distributed 
Ti with finite variance. 

Theorem 1: Let p{t) be as defined in equation (8) and Ti ~ Af(0, wifh > 0 a consfanf and 
^ < B < oo fox dll i, B & consfanf. Also, lef Ki be defined as in Secfion II-B and be independenf from 
Ti for all i. Then, limTv^oo ^Ar(f) = has fhe properties 

• ^00(213) — o> 

• Aoo{t) > 0 for t G (to — T,To), and Aoo{t) < 0 for t G (tq,tq + r) for some r < Tnz- 

• Aoo{t) is odd around t = tq, i.e. ^ 00(^0 + 0 = —^oo('ro — 0 for C ^ 0 

• Aoo(f) is continuous. A 

The properties ouflined in Theorem 1 will be key fo fhe synchronization mechanism we describe. The 
specific value of will be determined by our choice of fhe pulse-connection funcfion. Before we prove 
Theorem 1 in Secfion III-C.2 we develop and mofivafe a few importanf relafed lemmas. 
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1) Polarity and Continuity of Aoo{t): At time t = ri / tq, we have that 



1=1 


i=l 


where Mj(ri) = AmaxKip{Ti — tq — Ti). We have the mean of Mj(ri) being 


E{Mi{Ti)) = AmaxE{Ki) / p(ri - To - f)fTM)df, 


/ 


( 11 ) 


where /r. (i/)) is the Gaussian pdf 



(V’) 


2 


It is clear that the Mj(ri)’s, for different i’s, do not have the same mean and do not have the same 
variance since the two quantities depend on the a* value. Since the afs are characterized hy fa{s) 
(defined in Section III-A), we write the Gaussian distribution for T as 



and Mi{Ti) is in fact a function of s as well, denoted Mj(ri,s). Using fT{'f,s) and Mj(ri,s), the 
notation makes it clear that we can average over the a^’s that are characterized by /^(s). We use the 
results of Lemmas 1 and 2 to prove the polarity result for Aoo{t) in Section III-C.2. 

Lemma 1: Given the sequence of independent random variables Mi(Ti) with ti < tq, (ti)) = pi, 
and Var(Mj(Ti)) = af. Then, for all i. 


72 > /ii > 7i > 0 


( 12 ) 


a‘f < 73 < oo, 


(13) 


for some constants 71 , 72 , and 73 and 



i=l 


almost surely, where 




pin - To- 'lp)fTi'lp, s)d'ijjfa{s)ds. A 
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Lemma 2: Given the sequence of independent random variables Mj(Ti) with ti > tq, E{Mi{Ti)) = /r*, 
and Var(Mj(Ti)) = af. Then, for all i, 

72 < /ii < 7 i < 0 (14) 

a‘f < 73 < cx), (15) 


for some constants 71 , 72 , and 73 and 




N 


N^oo N 


2 = 1 


almost surely, where 


rcx-up 

r/(ri) = / E{Mi{n,s))fa{s)ds. A 

j ainu, 


The results of Lemma 1 and Lemma 2 are intuitive since given that p{t) is odd and the Gaussian noise 
distribution is symmetric, it makes sense for ^ 00 (f) to have properties similar to an odd waveform. Since 
the proofs of the two lemmas are very similar, we only prove Lemma 1. The proof can be found in the 
appendix. 

Knowing only the polarity of Aoo{t) is not entirely satisfying since we would also expect that the 
limiting waveform be continuous. The proof of Lemma 3 is once again left for the appendix. 

Lemma 3: Using p{t) in ( 8 ), 

1 ^ 1 ^ - 
^oo(f) = lim AmaxKip{t - To - Ti) = lim = rj{t) 

N^OO iV N^OO I\ 

2=1 2=1 

is a continuous function of t, where 


r]{t) 


/ E{Mi{t,s))fa{s)ds 

JoLion, 

rcxup roc 

AraaxE{Ki) / / p{t - Tq - 1 p)fT{'lp, s)d'll;fa{s)ds. 

j Oi-louj J—00 


A 
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2) Proof of Theorem 1: We can proceed in a straightforward manner to show that ^cx)('ro) = 0. For 


t = To, 

^ A K 1 ^ 1 

^n{to) = ^ ^ - Tq - Ti) = —J2AmaxKip{-Ti) = 

i=l i=l i=l 

where Mj = -AmaxKip{Ti). 

Since our goal is to apply some form of the strong law of large numbers, we first examine the mean 
of Mi. We have that E{Mi) = — AmaxE{Ki)E{p{Ti)). Furthermore, 

/ C50 

-oo 

since p{'f) is odd and fr^f’) is even because it is zero-mean Gaussian. Thus, E{Mi) = 0. 

We next consider the variance of Mp 


Var(M,) = E{Mf)-E\M,)=Al,,E{Ky{Ti)) 
= Al^axE{Kf)E{p^{Ti)) < < cx), 


where we have used the fact that E{Kf) < 1 and \p{t)\ < 1. 

From the preceding discussion we see that the Mfs are a sequence of zero mean, finite (but possibly 
different) variance random variables. From Stark and Woods [38], we know that if Var(Mj)/f^ < oo, 
then we have strong convergence of the Mfs: 


1 


N 


1=1 


Ah ^ E{Ah), 


with probability-1 as —> oo. But it is easy to see that 


Var(Mi) ^ 

2-^ j2 2-^ i 

i=l i=l 




max _ a '2 ^ 

^max g ^ ^^5 


SO the condition is satisfied. As a result, 


An{to) = ^ 0, 

i=l 

as N ^ oo. 

We have that Aoo{t) is continuous from Lemma 3. Thus, next we need to show that Aoo{t) > 0 
for f e (to — r, To), and Aoo{t) < 0 for f G (To,ro -|- r) for some r < r„^. We show the case for 
t = Ti £ {tq — r, To) by simply applying Lemma 1. Since Lemma 1 holds for all ri < tq, there clearly 
exists a r such Aoo(f) > 0 for f G (tq — t, tq). The case for t G (To,ro -|- r) comes similarly from 
Lemma 2. 
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Lastly, it remains to be shown that ^oo(i) is odd around t = tq. This, however, is evident from the 
form of T/(f). Since /r('0) s) is even in ijj about 0 and p(V’) is odd about 0, it is clear that p(t — To — 
'ip)fT{'>p,s)d'ip as a function of t is odd about tq. Thus, 7]{t) is odd around tq . This then completes the 
proof for Theorem 1. A 


IV. Asymptotic Time Synchronization 
A. The Use of Estimators in Time Synchronization 

In this work we want to show that as we let —> oo then we can recover deterministic parameters 
that allow for time synchronization. Such a result would provide rigorous theoretical support for a new 
trade-off between network density and synchronization performance. To simplify the study, we focus 
on the steady-state time synchronization properties of asymptotically dense networks. In particular, we 
develop a cooperative technique that constructs a sequence of equispaced zero-crossings seen by all nodes 
which allows the network to maintain time synchronization indefinitely given that the nodes start with 
a collection of equispaced zero-crossings. Starting with a few equispaced zero-crossings allows us to 
avoid the complexities of starting up the synchronization process but still allows us to show that spatial 
averaging can be used to average out timing errors. If we are able to maintain indefinitely a sequence 
of equispaced zero-crossing using cooperative time synchronization, then it means that spatial averaging 
can average out all uncertainties in the system as we let node density grow unbounded. This recovery of 
deterministic parameters is our desired result. Here, we overview the estimators needed for cooperative 
time synchronization. 

Let be the time, with respect to clock c^, that the Ah node sees its nth pulse. In dealing with the 
steady-state properties, we start by assuming that each node i in the network has observed a sequence 
of m pulse arrival times, that occur at integer values of t, m is an integer. Recall 

that 1 ..., ^ ■ is defined as a set of m pulse arrival times in the time scale of &. Therefore, 

even though i, ■ ■ ■, tn-m i occur at integer values of t (the time scale of ci), these values are not 
necessarily integers since they are in the time scale of Cj. Note also that in our model the pulse arrival 
time is a zero-crossing location. Using these m pulse arrival times, each node i has two distinct, yet 
closely related tasks. The first task is time synchronization. To achieve time synchronization, node i 
wants to use these m pulse arrival times to make an estimate of when the next zero-crossing will occur. 
If it can estimate this next zero-crossing time, then it can effectively estimate the next integer value of 
t. This estimator can then be extended to estimate arbitrary times in the future which gives node i the 
ability to synchronize to node 1. The second task is that node i needs to transmit a pulse so that the 
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sum of all pulses from the N nodes in the network will create an aggregate waveform that, in the limit 
as —> oo, will give a zero-crossing at the next integer value of t. This second task is very significant 
because if the aggregate waveform gives the exact location of the next integer value of t, then each 
node i in the network can use this new zero-crossing along with j, • • •, i to form a set of 

m zero-crossing locations. This new set can then he used to predict the next zero-crossing location as 
well as node Ts next pulse transmission time. Recall that determining the pulse transmission time is 
the joh of the pulse-connection function With such a setup, synchronization would he maintained 
indefinitely. The zero-crossings that always occur at integer values of t would provide node i a sequence 
of synchronization events and also illustrate how cooperation is averaging out all random errors. 

The waveform properties detailed in Theorem 1 play a central role in accomplishing the nodes’ task 
of cooperatively generating an aggregate waveform with a zero-crossing at the next integer value of t. 
From (10), if the arrival time of any pulse at a node j is a random variable of the form r^ + Ti, where tq 
is the next integer value of t and Tj is zero-mean Gaussian (or in general any symmetric random variable 
with zero-mean and finite variance), then Theorem 1 tells us that the aggregate waveform will make a 
zero-crossing at the next integer value of t. This idea is illustrated in Fig. 6. 



Fig. 6. Theorem 1 is key in explaining the intuition first illustrated in Fig. 1. The pulse p{t) is shown on the left 
figure, with tq = 1 and Amax = 1- On the right we have a realization of A]si{t) (N = 400), and we assume that 
Kj^i = 1 (no path loss) and Ti ~ Af(0, 0.01) for all i. As expected from Theorem 1, we notice that the zero-crossing 
of the simulated waveform is almost exactly at t = 1. 

Thus, for achieving time synchronization in an asymptotically dense network we need to address two 
issues. First, we need to develop an estimator for the next integer value of t given a sequence of m pulse 
arrival times that occur at integer values of t. We will call this the time synchronization estimator and let 
us write as the time synchronization estimator that determines the time, in the time scale of Cj, when 
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node i predicts it will see its nth zero-crossing. Two, we need to develop the pulse-connection function 
such that node i’s transmitted pulse will arrive at a node j with the random properties descrihed in 
Theorem 1. 


B. Time Synchronization Estimator Performance Measure 

Here we establish the conditions for estimating the next pulse arrival time, or equivalently the next 
integer value of t, given m pulse arrival times. These conditions apply most directly to the time syn¬ 
chronization estimator since we want to synchronize in some desired manner. The problem of 
synchronization is the challenge of having the ith node accurately and precisely predict when the next 
integer value of t will occur. In our setup, the reception of a pulse by node i tells it of such an event. 

Let us explicitly model the time at an integer value of t in terms of the clock of node i. Assume tq 
is an integer value of t and at this time, node i will observe its reth pulse. Thus, from (1) we have that 


= “*('^0 - Aj)-h T'i(ro). (16) 

The equation makes use of the clock model of node i (1) to tell us the time at clock Cj when node 1 is 
at To, where tq is an integer in the time scale of ci. We are also starting with the assumption that the 
zero-crossing that occurs at an integer value of t is observed by node i at this time. 

From (16) we see that the pulse receive time at node i, is a Gaussian random variable whose mean 
is parameterized by the unknown vector = [oj, tq, Aj]. Thus, to achieve synchronization node i will try 
to estimate the random variable using a series of m pulse receive times as observations (recall that 
m is known). Note that the observations are also random variables with distributions parameterized 
by -d. We want the time synchronization estimator of node i to make an estimate of denoted 
■ ■ ■ ^t'^n-rn,i) which is a function of past observations *at 

meets the following criteria: 




(17) 


argmin£c,^F;^ [(f^-..., 


(18) 


for all -d. The subscript r) means that the expectation is taken over the distributions involved given any 
possible 'd. The first condition comes from the fact that given a finite m, it is reasonable to want the 
expected value of the estimate to be the expected value of the random variable being estimated for all 
■d. As in the justification for unbiased estimators, this condition eliminates unreasonable estimators so 
that the chosen estimator will perform well, on average, for all values of r? [33]. The second condition is 
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the result of seeking to minimize the mean squared error between the estimate and the random variable 
being estimated for all 


C. Time Synchronization Estimator 

For the time synchronization estimator, node i will seek to estimate given v 

From (16), we see that T = j,..., is a jointly Gaussian random vector parameterized 

by r?. Recall that we assume 'l'j(f) is a zero mean Gaussian process with independent and identically 
distributed samples 'l'j(f) ~ AA(0,cr^), for any t. Also, since we’re assuming that the zero-crossings at 
node i occur at consecutive integer values of t, the random variable j is Gaussian with j ~ 
M{ai{To — m — Ai),a‘^) for some ?? = [o;*, tq — m, A,]. We also notice that 


= “*('^0 - m + 1 - Aj) = ai{To - m - Aj) + at. 


Since each noise sample is independent, we see that the distribution of T parameterized by ?? can be 
written as T AA(M, S) where 


M = 


at (to -m - Ai) 
aj(ro - m - Ai) + ai 
ai{TQ -m - Ai) + 2ai 


ai(To -m - Ai) + {m- l)ai 


and S = cj^I. 

As a result, for any m consecutive observations, we can simplify notation by using the model 


Y = H6» +W, 


(19) 


where Y = [yi .. .Y^Y = [f^_. 

e = 

with 




...G 


and 


n—m,i n—mH-1,2 * * * n—l,ii 

aj(ro -m- Ai) 
ai 


- 1 

_ 1 

1 

- 1 

1 - 

1 _ 

1 

1 - 


H 


-\ T 


111 ... 1 
0 1 2 ... m — 1 


and W = [VFi... WmY- Since 'l'i(f) is a Gaussian noise process, W ~ AA(0, S) with S = cr^I. 
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Using the simplified notation in (19), we want to estimate Ym+i, where Ym+i is jointly distributed 
with Y as 


Y 


M 

5 

S 0 

Ym+1 


9i + m02 


0 cj2 


Using this notation, we can rewrite the synchronization criteria as: 

Ee [Ym+i {Yi,Y2,...,Ym)]=Ee {Y^+i ) (20) 

argmiriy^^^Ee [(Yyn+i{Yi,Y2,Ym) - Y^+i)^], (21) 

where Ym+i is the estimator for Ym+i- 

Condition (20) implies that our estimate must he unbiased. Condition (21) is equivalent to 

avgmmy^^^Eg [{Ym+i{Yi,Y2, ..., Ym) - {9i + m6'2))^]. 

To see this equivalence, note that 

Eg [{Ym+l{Yi,Y 2 , ...,Ym)- Y^+l)'] 

= Eg [{Ym+l{Yi,Y 2 , ...,Ym)-{ei+ - Wm+l)^] 

= Eg[{Ym+i{Y^,Y2,...,Ym)-{ei + me2))^]+E[W^y^], (22) 

where the last inequality follows from the independence of Wm+i from all other noise samples. Since 
the distribution of of Wm+i is independent of 9, 

argminy ,^^^Eg |^(Tj, 2 _|_i (Y^, Y 2 , • • •, Ym) Ym-\-i) ] 

= argmin^^^^Eg [{Ym+i{Yi,Y2,... ,Ym) - {9i + m^s))^]. 

With these two conditions, from [33] we see that the desired estimate for Ym+i will be the uniformly 
minimum variance unbiased (UMVU) estimator for Eg{Ym+i) = 9i+ m02- 

Using the above linear model, from [23] we know the maximum likelihood (ML) estimate of 9, 9ml, 
is given by 

Oml = = (H^H)"^H'^Y. (23) 

This estimate achieves the Cramer Rao lower bound, hence is efficient. The Fisher information matrix is 
I{9) = ^^ 2 ^ and 9ml ~ AA(0, (T^(H^H)“^). This means that 9ml is UMVU. 

Again from [23], the invariance of the ML estimate tells us that the ML estimate for (j) = g{9) = 
9i + m02 is ^ml = 9iml + fn92ML- First, it is clear that ^ml = C9ml, where C = [1 m]. As a 
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result, we first see that EQ{(j)ML) = CEo{9ml) = + rn 62 so (pML is unbiased. Next, to see that (pML 

is also minimum variance we compare its variance to the lower bound. 



The extension of the Cramer Rao lower bound in [23] to a function of parameters tells us that 


Ee{\\g-g{e)f)>G{e)l-\d)G'^{e) 

with G{0) = {Vgg{ 6 ))'^. In this case, G{9) = [1 m] so the lower bound to the mean squared error is 



As a result, we see that 4>ml is UMVU. Since 0 ml is the desired estimate of where the next pulse 
arrival time will be, it is the time synchronization estimator. Thus, 


= C(H^H)-iH^Y. 


(24) 


Note that 



(25) 


has a variance that goes to zero as m —> oo. 

D. Time Synchronization with No Propagation Delay 

We now need to develop the pulse-connection function so that the conditions for Tj in Theorem 1 are 
satisfied. Recall we are developing the synchronization technique under the assumption of no propagation 
delay, i.e. 5{d) = 0. Given a sequence of m pulse arrival times, the time synchronization estimator V^\ 
given in (24) gives each node the ability to predict the next integer value of t. What remains to be 
considered is the second part of the synchronization process: developing a pulse-connection function 
X‘^- such that the aggregate waveform seen by a node j will have the properties described in Theorem 1 . 
Let us first consider the distribution of From (25), we have that 



Using (1), we can translate V^\(Y) into the time scale of ci as 

Kim = Mvzm) - A.) + 


which gives 
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This means that 


OY)~A'(ro,^(l + 


2(2m + 1) 


(26) 


m{m — 1) 

Under our assumption of 6 {d) = 0, any transmission hy node i will he instantaneously seen hy any node 
y. As a result, the random variable will he seen as the pulse arrival time at node j, in the time 

scale of Cl- 

Due to the assumption of no propagation delay, defining = V^\{Y) will give us the desired 

properties in the aggregate waveform. To see this, let us compare the distribution of X^^-(Y) to the 
assumptions of Theorem 1. Since tq is the ideal crossing time in the time scale of ci, we have 


K]i(y) = ro + T,. 

Therefore, we see that 

2 

Var(ri) = ^ 

where from Theorem 1 is 

\ m[m — 1) / 

We have shown that using the pulse connection function Ar^^-(Y) = V^\{Y) satisfies the conditions of 
Theorem 1. Thus, all the results of the theorem apply. 

As a result, we have established a time synchronization estimator V^\{Y) and a pulse-connection 
function X^-(Y). In the case of 6 {d) = 0, we have that X'^^{Y) = V^\(Y), or in the time scale of q, 
^ Ki i(^)- When each node in the network uses the pulse-connection function Y^v(Y) we have 
a resulting aggregate waveform that has a zero-crossing at the next integer value of f as —> oo. This 
fact follows from applying Theorem 1. Thus, we have an asymptotic steady-state time synchronization 
method that can maintain a sequence of equispaced zero-crossings occurring at integer values of t. An 
interesting feature of this synchronization technique is that no node needs to know any information about 
its location or its surrounding neighbors. 

1) Cooperation without Simultaneous Transmission and Reception: Before ending this section, let us 
comment on the assumption of simultaneous transmission and reception. One way to relax this assumption 
is to divide the network into two disjoint sets of nodes, say the odd numbered nodes and the even numbered 
nodes, where each set is still uniformly distributed over the area. Then, the odd nodes and the even nodes 
will take turns transmitting and receiving. For example, the odd numbered nodes can transmit pulses at 
odd values of t and the even numbered nodes will listen. The even numbered nodes will then transmit 
pulses at the even values of t and the odd numbered nodes will listen. With such a scheme, nodes do 


1 + 


Z[/m + ij 
mtm — 1 i 


cr 


v2 ’ 


(27) 
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not transmit and receive pulses simultaneously, but can still take advantage of spatial averaging. The odd 
numbered nodes will see an aggregate waveform generated by a subset of the even numbered nodes and 
the even numbered nodes will receive a waveform cooperatively generated by the odd numbered nodes. 
Let us take a more detailed look at this scheme. 

Aggregate signals generated by odd numbered nodes Aggregate signals generated by even numbered nodes 



Fig. 7. In the above figure, we assume tq is an even integer value of t and m = 3. Therefore, each even numbered 
node will turn on its receiver to receive the aggregate signal arriving at times tq — 5, tq — 3, and tq — 1. Using these 
three received times, it can then estimate the time of tq. Thus, the aggregate signal occurring at tq is cooperatively 
generated by the even numbered nodes and is received by the odd numbered nodes. 


In Fig. 7 we assume that tq is an even integer value of t and use m = 3. Each even numbered node 
will use the aggregate signals occurring at tq — 5, tq — 3, and tq — 1 to estimate tq and cooperatively the 
even nodes will generate the aggregate signal at tq. The odd numbered nodes will then use the aggregate 
signals occurring at tq — 4, tq — 2, and tq to generate the aggregate signal at tq + 1. Therefore, the 
odd and even numbered nodes can take turns transmitting and receiving signals and nodes never need to 
simultaneously transmit and receive. 

Of course, such a setup would require a modification of the estimators used by the nodes. Nodes 
will receive a vector of m observations Y with Y[/ + 1] = ai(ro + 1 — 2(m — 1) — Aj) + 'Fj for 

/ = 0,1,..., m — 1. With such a mechanism, the H matrix in equation (19) would change to 

r -\T 


H = 


1 1 1 
0 2 4 


1 

2(m — 1) 


and 6 becomes 


01 


ai{To + 1 - 2m - Aj) 

_ 02 _ 


0^2 


To estimate the location tq in the time scale of q, we can proceed as in Section IV-C: 

Oml = (H^S-^H)-^H^S-^Y = (H^H)-^H^Y 
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will be distributed 9ml ~ AA(0, and 9ml is UMVU. This leads to the UMVU estimate 

4>ml = C9ml, where C = [1 2m - 1], and E{^ml) = CE{9ml) = 9i + (2m - 1)02- In this case, 
the variance of (pML will be Var 5 )((^ML) = C(T^(H^H)"^C^, and thus we have that 


= ^ml ~ AA( ai{To + 1 - 2m - Ai) + (2m - l)ai, 


(T^(2m + l)(2m — 1) 
m(m — l)(m + 1) 


Converted to the time scale of ci we have 


C.(Y)~V 1 + 


CJ 


a: 


(2m + l)(2m — 1) 
m(m — l)(m + 1) 


(28) 


Comparing equations (26) and (28), we see that they have the same form. As a result, we can again set 
^niO^) ^ achieve cooperative time synchronization. 


V. Time Synchronization with Propagation Delay 

We now extend the ideas of cooperative time synchronization to the situation where signals suffer not 
only from pathloss but also propagation delay. It turns out that the effect of propagation delay can also 
be addressed using the concept we have been using throughout this paper — averaging out errors using 
the large number of nodes in the network. 

In this section, we use the pathloss and propagation delay model detailed in Section II-C. We introduce 
a time delay function 5{d). For generality, we explicitly model a multi-hop network where we have a 
K{d) function that is zero for d greater than some distance R, i.e. K{d) = 0 for d > R. Such a model 
implies that the aggregate signal seen at any node j is influenced only by the set of nodes inside a 
circle of radius R centered at node j. With this we can effectively divide the network into two disjoint 
sets, a set of interior nodes and a set of boundary nodes. An interior node j is defined to be a node 
whose distance from the nearest network boundary is greater than or equal to R. A boundary node is 
thus defined to be a node that is a distance less than R away from the nearest network boundary. 

We make this distinction since the synchronization technique for each set of nodes is different. Please 
note that if a pathloss function where K{d) = Q for d> R is unreasonable, then we simply choose R to 
be infinite and consider all nodes in the network to be boundary nodes. 

Using the propagation delay model, Dj^i will obviously modify the general received aggregate wave¬ 
form seen at any node j. In fact, equation (7) will now be written as 

N . „ 

AAW = E e9) 

i=l 

For N large, this model will give an accurate characterization of the aggregate waveform seen at node j. 
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A. Conceptual Motivation 

From equation (29), it is clear that the aggregate waveform will not have a zero-crossing at tq for every 
node j because of the presence of the Dj^i random variables. Therefore, to average out propagation delay, 
the idea we employ is to have each node introduce a random artificial time shift that counteracts the effect 
of the time delay random variable. More precisely, we want to introduce another random variable Dfi^ 
such that D+ Dj will have zero mean and a symmetric distribution. At the same time, we assume each 
node knows K{-) and 6{-) and will also introduce an artificial scaling factor Kfi^ = K{6~^{—Dfix)) 
to simplify the analysis of the aggregate waveform. This means that instead of using the scaling factor 
Ai = Amax/N, each node i will scale its transmitted pulse by A, = AmaxKjix/N. For the motivation 
in this section, let us assume that node j is an interior node. 

To find fhe distribution of Dfi^, we consider the following. Dj has density fojix) and let be 

the density of Djix- Since Dj and Djix are independent, we know that the density of Dt = Dfix-\-Djj, 
fDxix), will be the convolution of fojix) and /£)^.^(x). Therefore, by the properties of the convolution 
function, if we set fuft^^ix) = fDj{—x), then we have that fDxix) is symmetric, i.e. forix) = fDT{~x)- 
As well, since Dj has finite expectation, it is easy to see that E{Dt) = 0. 

Given a sequence of m zero-crossings that we know to be occurring at integers of t, we can still use 
V^\{Y) (from (24) in the time scale of node 1) as the time synchronization estimator. However, with 
propagation delay, the pulse-connection function will now be = V^\{Y)+Dfix = Xo+Ti+Dfix- 

With Dfix and Kfix included, we can rewrite equation (29) as 

= E 2=£:&A2£p(t _ - r, - - D,.,). (30) 

i=l 

It is important to see that since Dj has the same distribution for all interior nodes j, equation (30) holds 
for every node j that is an interior node. This means that for the network to cooperatively generate 
the waveform in (30) each transmit node i needs to have the following additional knowledge: (1) the 
distribution of Dfix whose density is fDfi^{x) = fn^i—x), where j is an interior node, and (2) the 
functions K{-) and (5(-) to generate Kfix- With this knowledge, we can use equation (30) to study the 
aggregate waveform seen at any interior node j. In fact, we find that the aggregate waveform has limiting 
properties that are similar to those outlined in Theorem 1. These properties are described in Theorem 2. 

Theorem 2: Let p{t) be as defined in equation (8) and Tj ~ AA(0, ^) with > 0 a constant and 
^ H < (X) for all z, H a constant. Kjf and Djf are defined as in Section II-C and Dfix with 
density /ofiAx) = fn^i-x) is independent from Djf. Kfix = K{6~^{-Dfix)) and let Djf, Dfix, 
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and Ti be mutually independent for all i. Then, for any interior node j with as defined in (30), 

limAr^oo^j"Ar(f) = has the properties 

• ^loci'^o) = 0 , 

• ^^' 00(0 is odd around t = tq, i.e. Af^^{To + 0 = - 0 for ? > 0. A 

The proof of Theorem 2 is left for the appendix. 

From the arguments so far, it seems that time synchronization with delay, at least for interior nodes, 
can be solved simply by modifying the pulse-connection function X^^^(Y) and changing the scaling 
factor to Ai = AmaxKfi^/N. Theorem 2 tells us that the limiting aggregate waveform makes a zero¬ 
crossing at the next integer value of t and the waveform is odd. Thus, we can use this zero-crossing as 
a synchronization event and maintain synchronization in a manner identical to the technique used in the 
situation without propagation delay. This, however, unfortunately is not the case. In order to implement 
the above concept, we need to find fhe random variable, Dj-^, in the time scale of a, that corresponds 
to Dfix such that 


(ICi (Y) + = 


v;-(Y) + - 5. 


a,; 


+ Aj 


= VUY) H- ^ 


This means that we need = Dji^. However, each node i cannot find that satisfies fhis 

since it does not know its Oj. 


B. Time Synchronization of Interior Nodes 

Since the fth node does not know its own value of Oj, to do time synchronization with propagation 
delay we can have each node estimate its Oj value. However, this estimate will not be perfect and we may 
no longer have the symmetric limiting aggregate waveform described by Theorem 2. This means that the 
center zero-crossing might occur some e away from tq, tq an integer value of t. However, steady-state 
time synchronization can be maintained if the network can use a sequence of m equispaced zero-crossings 
that occur at f = tq — m + e, tq — m -|- 1 -|- e, tq — m -|- 2 -|- e, ..., tq — 1 -|- e, where tq is an integer 
value of t, to cooperatively generate a limiting aggregate waveform that has a zero-crossing at tq -|- e. 
In such a situation, the network will be able to construct a sequence of equispaced zero-crossings and 
maintain the occurrence of these zero-crossings indefinitely. The idea is the same as in the case without 
propagation delay, but the only difference here would be that the zero-crossings do not occur at integer 
values of t. Let us give a more formal description of this idea. 
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Using notation from Section IV-C, we start with the assumption that each interior node i has a sequence 
of m observations that has the form 


aj(ro - m + ^ + e - Aj) + (31) 

where ( = 0,l,...,m — 1 and e is known. To develop the time synchronization estimator V^\{Y) and 
the pulse-connection function X^\(Y), we consider the observations made by each node. If we assume 
that each node knows the value of e, the vector of observations can be written as in (19) 


Y = H6» +W, 


where the matrix H in this case is 


H 


-\ T 


1 1 1 ... 1 
e 1 + e 2 + e ... m — 1 + e 
Using this model, we can follow the development in Section IV-C to find the the time synchronization 
estimator 


U„^UY,e) = C(H^^'H)~iH^Y, 




(32) 


where C = [1 m]. This estimator will give each node the ability to optimally estimate the next integer 
value of t. Note that the variance of the time synchronization estimator is 


Var0(U„=^(Y,e)) = Ca\WYi)-^C^ = a 


TxTX-lr^T 


2 / 2(2m + 1) 12e(e — 1 — m) 


m{m — 1) (m — l)m{m + 1) 

Using the time synchronization estimator, we can choose the pulse-connection function as 


(33) 


X;:v(Y) = V::,iY,e)+a,Dfi, = v::,(Y ,e) + D%, 


where each time node i makes the estimate V^^(Y,e) it also estimates dj as 


di = C(H^ H)-^H^Y, 


(34) 


C = [0 1]. We find that di ~ M{ai, 12a'^/{{m — l)m{m + 1))). Since, from Section V-A, we know 

we want D^fi^/ai = Djix, we have set = diDfix- Notice that since is simply a realization 
of Dfix multiplied by node i’s estimate of a*, node i can use the realization of Dfix and find Kfix = 
K{S-^{-Dfix)). 

With our choice of X^\{Y) in (34), we see that 

iV:}i^x)+D%r = V:i,iY,e) + Z,Dfix = ro + Ti + Z^Dfix, 
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where Zi ~ Af{l, 12a'^/{af{m — l)m{m + 1))), and tq + Tj = e). Because of the random factor 

Zi, we see that Dt = ZiDfi^ + Dj^i is no longer a symmetric distrihution. As a result, the limiting 
aggregate waveform 



may not have a zero-crossing at f = tq. 

Thus, if we can find an e such that each node i using a set of observations of the form (31) allows the 
network to cooperatively generate the waveform in (35) that has its zero-crossing occurring at f = T^+e (in 


the time scale of ci), then we have steady-state time synchronization. This is because the network would be 


able to use a sequence of m observations to generate the next observation that gives the same information 
as any of the previous observations. Thus, by always taking the m most recent observations, the process 
can continue forever and maintain synchronization. Each node i would need to know distribution of Dfi^, 
the value of e, and the functions K{-) and (5(-). Therefore, we find that steady-state time synchronization 
of the interior nodes is possible under certain conditions. As a note, no interior node needs to know any 
location information. 

C. Time Synchronization of Boundary Nodes 

Before we consider the synchronization of boundary nodes, we note that the key requirement for each 
boundary node i is to have a pulse-connection function given in equation (34). The reason that this must 
be the pulse-connection for every boundary node i is because the analysis for the interior nodes assumes 
that the aggregate waveform seen by any interior node j is created by pulse transmissions occurring at a 
time determined by (34). Since the aggregate waveform seen by some interior nodes are created by pulse 
transmissions from boundary nodes, each boundary node must have the appropriate pulse-connection 
function. This requirement, however, proves to be extremely problematic and reveals a limitation of the 
elegant technique of averaging out timing delay when we come to boundaries of the network. 

The problem comes because Dfix + Dj^i already does not have a symmetric distribution if j is a 
boundary node. Recall that = fDj{—x) when j is an interior node and fDj{x) = foiix) when 

j and I are both interior nodes. However, fi). (x) / /d, (x) when j is an interior node and I is a boundary 
node. As a result, Dfix + Djf is no longer symmetric if j is a boundary node. In fact, it is clear that the 
distribution of Dfix + Djf is a function of node j’s location near the boundary. Because of this additional 
asymmetry, let us assume for a moment that the sequence of zero-crossings observed by boundary node 
i occur Ci away from an integer value of t. That is, if every node in the network, including the boundary 
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nodes, transmitted a sequence of pulses where each pulse was sent according to (34), then boundary node 
i would observe the sequence of observations 

ai(ro - m + / + ej - Aj) + 4'j, (36) 


where / = 0,1,... , m — 1 and e* is known. 

This boundary node i could then use the time synchronization estimator given by (32) but where the 
matrix H is now replaced with Hj 


H, 


1 1 1 
Ci 1 + Ci 2 + Cj 


Thus, for this boundary node i we have 


1 


T 


m-l + ei 




(Y,ei) = C(HfH,)-iHfY, 


(37) 


In this case, however, the variance of the time synchronization estimator depends on ei 


\^xg{VZ{Y,ei)) = a^ 


2 / 2(2m + 1) 12ej(ei — 1 — m) 


+ 


m{m — 1) {m — l)m{m + 1) 

The fact that the variance depends on e* is the root of the problem. The pulse-connection function 


(38) 


^::t(Y) = ^(Y,e,) + a,Z7y,., 


(39) 


is not the same as that given by (34). 

To correct for this, we can make the strong assumption that each boundary node i knows is own a*. 
We address the reasoning behind this assumption in Section V-D. If we use this assumption, then each 
boundary node i can get an observation sequence of the form (31) simply by adding ai{e — e*) to each 
of the m observations of the form given in (36), where we assume that node i knows both e and e*. 
With such an observation sequence, boundary node i will have the time synchronization estimator (32) 
and, more importantly, the pulse-connection function (34). Thus, maintaining time synchronization for 
the case of propagation delay would be possible. 

What we have then is that boundary node synchronization would require only the boundary nodes 
to know their Oj parameters. With this strong assumption only for the boundary nodes, the network 
is effectively synchronized. Even though the boundary nodes do not see the same zero-crossing as the 
interior nodes, they can calculate this time and thus have all the required synchronization information. 
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D. The Boundary Node Assumption 

The assumption that each houndary node i knows ai is a strong assumption. Even though the fraction of 
nodes that are boundary nodes is small for multi-hop networks requiring many hops to send information 
across the network, we helieve that the assumption is still very artificial. There are two reasons that we 
make the assumption for the presentation of results on time synchronization with propagation delay. 

First, the assumption allows us to give an elegant presentation of the main concept of this paper which 
is to use high node density to average out errors in the network. Throughout this work we have used high 
node density to average out inherent errors present in the nodes. We were able to average out random 
timing jitter that is present in each node and provide the network with a sequence of zero-crossings 
that can serve as synchronization events. We then applied this technique to averaging out the errors 
introduced by time delay. To this end we were partially successful in that the interior nodes can average 
out these errors assuming the boundary nodes have additional information. But this is of interest since 
the goal of this paper is to understand the theory of spatial averaging for synchronization and discover 
its fundamental advantages and limitations. 

Second, the problem encountered at the boundaries is one that opens up an entirely new area of study 
which is the target of our future work. The issue that we encounter is that the waveform seen by some 
nodes in the network will have a zero-crossing that is shifted from the ideal location. This implies that 
different nodes will observe different zero-crossings. Furthermore, these zero-crossings will now evolve 
in time since we do not have the same observations over the entire network. This problem is similar 
to what we encounter if we consider finite sized networks. For finite N, the zero-crossing location will 
be random and thus introduce another source of error. As well, different nodes will see different zero¬ 
crossing locations. Therefore, we will turn our attention to the case of finite N and develop a different 
set of tools that will be needed to understand what types of synchronization are achievable under the 
situation where zero-crossing locations evolve in time. Using this understanding, we hope to return to 
the issue of propagation delay in asymptotically dense networks and characterize the behavior of the 
network. 


VI. Conclusions 

To conclude, we revisit the scalability issue under the light of work developed in this paper. 
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A. The Scalability Problem Revisited 


In the Introduction (Section I-B.2), we mentioned that most existing proposals for time synchronization 
suffer from an inherent scalability problem. The problem with those existing proposals lies in the fact that 
synchronization errors accumulate: if node 2 can synchronize to node 1 with some small error, and node 
3 can synchronize to node 2 with the same small error, these errors accumulate, and the synchronization 
of node 3 to node 1 is worse. Therefore, synchronization error increases with the number of hops in the 
network, and this problem is especially apparent in the regime of high densities. To make these ideas 
precise, we first determine the maximum number of hops over which synchronization information must 
travel and then study how the error in a generic pairwise synchronization mechanism depends on this 
number of hops. 

1) An Estimate of the Maximum Number of Hops: To obtain an estimate for the maximum number 
of hops in a network in the regime of high densities (fixed area, N —> oo), we approximafe fhe 
fransmission range of a node by fhe minimum required fransmission disfance, d^, to maintain a fully 
connected network with high probability. From [15], we have that for N nodes uniformly distributed 
over a [0,1] x [0,1] square, the graph is connected with probability-1 as —> oo if and only if each 
node’s transmission distance c^tv is such that 

2 log N + eN 

= -TT-, 


for some e^r ^ oo. Let us, therefore, approximate as 

11 log N 

Thus, — l/dM — and thus ^ oc as N —* oo. 

2) Synchronization Error Over Multiple Hops: Now, we assume there are nodes arranged in a 
linear ordering, numbered 1 to To synchronize, each node i forms an estimate of its own Oj, based 
on m pulses transmitted from node i — 1. As before, node 1 will have the reference clock ci(f) = t. 

Node 1 starts by sending m pulses at times ti + / for / = 0,1,... , m — 1. As a result, node 2 will 
get a vector of observations Y 2 , where Y 2 [l] = q; 2 (ti — ^ 2 ) + '['2 and the (/ + l)th element of Y 2 is 
Y 2 [( + 1] = (X 2 {ti — A 2 ) + la 2 + 'I' 2 - This is similar to the situation we had in (19) and we can therefore 
estimate 02 using 


02 = C(H^H)-1h^Y2, 

where C = [0 1]. We find fhat 0:2 ~ Jf{a 2 , 12cj^/((m — l)m(m + 1))). 
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Node 2 will now transmit m pulses at times, in terms of C 2 , T 2 + lci 2 , for / = 0,1,... , m — 1. Note 


that a 2 is now a fixed value since node 2 has estimated a 2 - In terms of ci, these pulses occur at 


(t 2 + la2Y" = 


T2 + la2 - 'I'2 
02 


+ Aa 


I 7 “2 

T2 + I - 

02 


02 ’ 


for f = 0,1,..., m — 1, where T 2 = (r 2 /o 2 ) + A 2 . Thus, if we translate these times into the time scale 
of C3, we will have the vector of observations, Y 3 , made hy node 3. We find fhaf fhe (Z + l)fh elemenf 
of Y 3 is 


Y3[Z + 1] = 03((r2 + I 


02 

02 


T'o 

— )- A3) + T'3 

O 2 


AA( 03(r2 - As) + Z03 —+ 


,an 


02 


Oio 



This vector of observations is of the form 


where 


with 


9 = 


H 


Y 3 = H0 + W, 

03(r2 - A3) 


h 

1 

- 1 

_ h _ 

1 

1 - 


1 1 1 
0 12 


03^ 

02 


1 

m — 1 


n T 


and W = [hFi... Wmf. W ~ M{0, S) with S = cj2(^ + l)l. 
With this vector of observations, we can use the estimator 


03 = C(H^H)“1h^Y3, 


where C = [0 1]. We find fhaf 

03 ~ 


AA(o3^r 


12 fj 2 


«2 ’ ((^ — l)m(m + 1 )) ''02 

If we continue fhis reasoning, we find fhaf 


(3+1) 




■M 






12 cj 2 


■( 


o 


tN 


1) 


(40) 


- l)m(m + 1 )) 

will be fhe esfimafe of node ij^. 

From fhe above analysis, we see fhaf each node Ts esfimate suffers from jitter variance of the same 
form. However, there is an accumulation of error because node i’s estimate has a mean that is dependent 
on node Z—I’s estimate. As a result, if node Z—1 has some small error, then that error will propagate to the 
estimate of node Z. A good way to see this is if we consider the special case where a 2 = = ■ ■ ■ = 1. 
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This is the case where the clock frequencies are the same, hut nodes do not know this. In this case, we 
find that node f^r’s estimate can he written as 

t-N 

= 02 + ^ Wi, > 2 

i =3 

where Wi ~ AA(0, 24cj^/((m — l)m(m + 1))). This is intuitively obvious because node z’s estimate cij 
will be the mean of the Gaussian random variable dj+i. Therefore, it is obvious that the error variance 
grows linearly with the number of hops. In fact, this behavior is observed in experimental work. With 
Reference Broadcast Synchronization (RBS), from [8] the authors find fhaf fhe synchronizafion error 
variance of an hop pafh is approximafely where is fhe one hop error variance. Therefore, 

we have fhaf fhe synchronizafion error befween our fwo nodes will grow linearly as which 

is sfricfly monofonically increasing. As a resulf, as ^ oo, we have fhaf synchronizafion error will 
grow unbounded. 

This scalabilify problem, however, can pofenfially be avoided using cooperative time synchronization 
as —> oo. This is because in fhe limif of infinife densify, fhe cooperafive lime synchronization technique 
allows every node in fhe nelwork lo see a sel of identical equispaced zero-crossings. As a resulf, in sfeady- 
sfale fhe synchronizafion error does nol grow across fhe nefwork. This comes abouf by using fhe high 
node densify lo average oul random liming errors. Thus, we find fhaf cooperative lime synchronization 
has very favorable scalabilify properties in fhe limif as A^ —> oo. 

B. Network Density and Synchronization Performance Trade-Off 

The cooperafive synchronizafion technique described in fhis paper provides us deferminisfic parameters 
fhaf we can use for lime synchronizafion in fhe limif as node densify grows unbounded. In facl, as fhe 
node density grows, fhe observalions fhaf can be used for synchronizafion improve. This means fhaf 
our cooperafive synchronizafion technique provides an effective Irade-off befween nelwork density and 
synchronizafion performance. Such a Irade-off has nol exisled before and will provide nelwork designers 
an addifional dimension over which fo improve nelwork synchronizafion performance. 

The fundamenlal idea behind cooperafive lime synchronizafion is fhaf by using spalial averaging, fhe 
errors inherenl in each node can be averaged oul. By using observations fhaf are an “average” of fhe 
information from a large number of surrounding nodes, synchronizafion performance can be improved 
due fo fhe higher quality observations. 

From fhis poinl of view, if is clear fhaf fhe parlicular technique described in fhis paper is bul one 
example of using spalial averaging fo improve synchronization. Olher techniques can also be developed 
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using spatial averaging. For example, nodes may not necessarily have to send odd-shaped pulses and 
use zero-crossing observations. Even though this setup takes advantage of the superposition of pulses, it 
has its drawbacks. To keep the signals in phase, the jitter variance will limit the maximum frequency at 
which signals can be sent. Instead, nodes may transmit ultra wideband pulses. If the nodes surrounding 
a particular node j each transmit an impulse at their estimate of an integer value of t, then due to timing 
errors in the surrounding nodes, node j will see a cluster of pulse arrivals around this integer value of 
t. Node j can then take the sample mean of this cluster of pulses and use that as an observation, just 
like we used the zero-crossing as an observation in this paper. This idea is illustrated in Fig. 8. Such 
a technique based on ultra wideband pulses will also provide similar scalability properties. As a result, 
cooperative time synchronization really describes a class of techniques that can take advantage of spatial 
averaging to improve synchronization performance. 


Pulse cluster 


time 



Aggregate waveform 


Fig. 8. Clusters of ultra wideband pulses can be used for cooperative time synchronization. In the the top figure, 
we illustrate the clusters of pulses around integer values of t. As the number of nodes increase, the sample mean 
will converge to the integer value of the reference time. This idea is parallel to the use of zero-crossings shown in 
the bottom figure. 


C. Future Work 

With the goal of developing practical cooperative synchronization mechanisms, two keys areas of 
interest are cooperative synchronization in finite-sized networks and algorithm development. First, the 
analysis of performance for finite-sized networks is very important. Determining when the asymptotic 
properties presented in this work are good predictors of performance in networks that may be large but 
still finite in size is important in terms of bridging the gap between our proposed ideas and practical 
systems. Preliminary, simulation-based work along these lines can be found in [19]. Second, developing 
practical techniques for cooperative time synchronization is essential for implementing spatial averaging 
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in real networks. Along these lines, one area of interest is determining what types of pulses should he 
used, i.e. odd-shaped pulses or ultra wideband pulses. 

Furthermore, the ideas in this paper suggest a few other areas of interest for future work. One is the 
issue of distributed modulation methods. If we have the ability to generate an aggregate waveform with 
equispaced zero-crossings, by controlling the location of these crossings we can modulate information 
onto this waveform and use it to communicate with a far receiver. Preliminary work along these lines 
can be found in [20]. Another issue is to study how the idea of spatial averaging that is so prevalent in 
this work contributes to synchronization that is observed in nature. 
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Appendix 

Proof of Lemma 1, To show (12), we consider 


E{AmaxKip{Ti - To -Ti)) = AmaxE{Ki)E{p{Ti - Tq - Ti)) 

= AmaxE{Ki) j p{ti - To- 
= -AmaxE{Ki) I p{ii - (n - To))fTiWdlp 

Since ri < tq, we have that ri — tq < 0 implying that p{'ip) is shifted to the left and the zero-crossing 
of pi'ip) occurs at a negative value. p{'ip) is odd about its zero-crossing and fTi{'>P) is symmetric about 
zero and strictly monotonically increasing on (—oo,0] for all positive finite variance values. Thus, it is 
clear that f p{\j) - (n - To))fTiWdip < 0 which makes E{AmaxKip{Ti - tq - Tj)) > 0 . 

Now, the expectation will vary with the variance of Ti and the variance will range from a positive 
upper bound of < B to a positive lower bound of where recall that is a value 

determined by our choice of the pulse connection function. If we consider f p{'il; — (ri — To))fTi {'4’)d'ip to 
be a function of the variance of Ti, then we see that it is bounded and continuous on the compact domain 
Since we showed in the previous paragraph that E{AmaxKip{Ti — To — Ti)) > 0 
whenever Ti has a nonzero finite variance, clearly E{AmaxKip{Ti — To — Ti)) > 0 when Var(rj) G 
Thus, it is clear that 71 and 72 exist and (12) is shown. 

To show (13), we consider 


\aT{AmaxKip{Ti - To - Ti)) = E{Af^^^Kfp^{Ti - To- Ti)) - E^{AmaxKip{Ti - To- Ti)) 


< Ai,,,E{Kt)E{p\T^-T0-T,)) 


< Al,,E{Kl) 


< A: 


— max 


where the second to last inequality follows from the fact that E{p‘^{ti — tq — Ti)) is upper bounded by 
1. The last inequality follows since E{Kf) < 1 by the fact that 0 < ifj < 1. Thus, we have shown (13). 

Next we define Sn = Mi(ri) -f • • • -f M„(ri) and rUn = E{Sn) = pi + ... + Pn. From [10] we have 
the following theorem 

Theorem 3: The convergence of the series 
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implies that the strong law of large numbers will apply to the sequence of independent random variahles 
Mj(Ti). That is, again from [10], for every pair e > 0, (5 > 0, there corresponds an N such that 

\Sn - rrir, 


Pr 


n 


<e; n = A^, A^ + l,...,A^ + r>>l —5 


for all r > 0 . 


A 


We have shown (13) so we have cr? < 73 < oo. Thus 


N 2 


N 


lim V 4 < lim 

^ 7 ,^ N^OC ^ I 

2=1 


N^oo I 
2 = 1 


73 

2 = 73 


TT 

Y 


and we have convergence hy the direct comparison test. Therefore, we can apply Theorem 3 and get that 
for any pair e > 0, 5 > 0, we can find an N such that 


Pr 


n 


m, 

n 


<e; n = A, A^ + l,...,A^ + r>>l —5 


(41) 


for all r > 0 . 

By (12) we have that 72 > /Tj > 71 > 0. Thus, we can clearly see that 

TTlr, 


n 


> 71- 


Furthermore, since we keep the function fa{s) constant as we increase the number of nodes in the 
network we get that mnjn converges to a constant rj{Ti) given by 


7]{ti) = A 

max EiKi) 


/ OO 

pin 

-OO 


-To- s)d'il}fa{s)ds 


f 

•J Oil 


E{Mi{Ti,s))fa{s)ds. 


The above expression comes from the fact that since each = i?(Mj(ri)) is a function of a*, mn/n will 
converge to the average of the jii over fa{s), the function that characterizes the set of ctj’s. Therefore, 
given any e, we can find an N' such fhat 

rUr, 


n 


- 7(n) 


< e 


(42) 


for all n > . Note that since (mn/n) > 71 , we have that ^(ti) > 71 . Since 

Sn '^n 
n n 


— -7(n) 
n 


< 


+ 


rrin , , 

- 

n 


using (41) and (42) we have 
Pr 


— -7(n) 
n 


< 2e; n = iV", N'^ + 1,..., N'^ + r } > 1 - 6. 
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for all r > 0, where N^' 


max{iV, A^'}. Thus, we have 


I™ ^ > 0 


N 


N—too N 


i=l 

almost surely. This completes the proof of Lemma 1 . A 

Proof of Lemma 3. First, we start hy finding an analytical expression for |Aoo(f) — ^cx)(fo)|- From 
the proof of Lemma 1 we have that 


A. 


ra^p noo 

(t) = AmaxE{Ki) / / p{t-To- s)dil:fa{s)ds 

^ OLlow —00 


Therefore, |Aoo(f) — ^ 00 (A)! can he written as 


l^oo(f) ^C»(fo)| 

rct^p poo 

= \AmaxE{Ki) / / \p{t - To - - p{to -To- V')]/r('i/’, s) fa{s)dll;ds\ 

^ Oilow —00 
pOLup poo 

< Amax / / \p{t -To-i)) - p{to -To- '0)|/t(V’, s) fa{s)d'lpds 

J OLiovj J— 00 

pa^p pTm+t„-To+\t-t„\ 

= Amax / \p{t - To - Ip) - p{to - To -'lp)\fT{i’,s)fa{s)d'll;ds, 

J Qlo-w d —Tnz+to-To — \t—ta\ 

where E{Ki) < 1. The change in the limits of integration in the last equality comes from the fact that 
p{i -To-lp) - p{to - To-pj) =0 outside of 'P € [-Tm + to - To - \t - to\,Tnz + to - Tq + \t - to|]. 
This is the maximum interval over which p{t — To — ip) — p{to — To — ip) can he non-zero. There is no 
need to take the absolute value of fT{ip,s) and fa{s) since they are always non-negative. 

Our second step is to hound the inner integral. Before doing so, we first show that the inside integral is 
in fact Riemann integrahle. For any given t and to, the inside integral is taken over a closed interval. Over 
a closed interval, we know from Strichartz [39] that any hounded function that is continuous except at a 
finite number of points is Riemann integrahle. Furthermore, also from [39] we know that the sums and 
products of continuous functions are continuous. As well, if a function is continuous then the absolute 
value of that function is also continuous. p{t) has at most D = 3 locations at which it is discontinuous 
and over any open interval not containing a discontinuity, p{t) is uniformly continuous since q{t) is 
uniformly continuous. fT{ip,s) has D' = 0 discontinuities in ip for an given s since it is Gaussian for 
any s. And since s G [o:iow,ctup], \fT{ip,s)\ < Gt for all ip and s {Gt occurring when ip = 0 and 
s = ttup)- Thus, since p{t) and friip, s) are continuous except at a finite number of points, we see that 
for given s, t, and to 

\p{t -To -Ip)- p{to -To- 1p)\fT{lp, S) 
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is also continuous in except at a finite number of points (at most D' + 2D points). This function is also 
bounded since the product of two bounded functions is bounded. As a result, we see that the integral is 
Riemann integrable over any closed interval. 

We now proceed to bound from above the value of this integral by first bounding the maximum value 
of the integral assuming no discontinuities and then introducing another term that bounds the maximum 
area contributed by the discontinuities. If we ignore the discontinuities and assume p{t) is uniformly 
continuous, for any mi > 0 there exists a n > 0 such that 


\t - fo| < - 
n 


\p{'t) -p{to)\ < —, 

mi 


for all t and to- As a result, p{t — To — ij)— p{to — To — i/j) can be made as small as desired by choosing 
the proper n thus giving us p(t — To — V’) “ p(^o — Tq — V’) < I/ttii for all ijj for an appropriate choice 
of n. 

Furthermore, we note that |p(f — Tq —'i/’)|/T('0, s) < Gt because \p{t)\ < 1 and \fT{'il^,s)\ < Gt- The 
maximum possible jump at a discontinuity in the function \p{t — To — ip) —p{to — Tq — '0 )|/r('!/’) s) is thus 
2Gt and for any |t —to|. the maximum area contributed by each discontinuity is 2GT\t — to\- As a result, 
for all D'+2D discontinuities, the maximum area contribution will be no more than 2GT\t—to\{D'+2D). 

We can, therefore, bound the inner integral as 


/ 


Tnz+to—ro + \t — to\ 


-T„z+tc-To-\t-to\ 
< 


\p{t -To-'tp)- p{to -To- S)dll) 


T,^z+to-To + \t-to\ n 

—d^p + 2GT\t - to\{D' + 2D) 

— Tnz+to-T 0 — \t — to\ 


/ 

= - i‘2'Tnz + 2|f — to\) + 2GT\t — to\{D' + 2D) 

mi 

= 2- Tnz + 2- \t — to\ + 2GT\t — to\{D' + 2D), 

mi mi 


where |f — to\ < Ijn. 

What we have is that if |f — to| < l/u- then 


l^oo(f) 2 loo(fo)| 

< ^max [ [2 - Tnz + 2 - \t — to\ + 2GT\t — to\{D' + 2D)'\ fot{s)ds 

Jaiou, \ mi mi J 


< A 


^ -^max 


Ga{oiup — Oilow) (^2 - Tnz + 2 - |f — fo| + 2G7’|f — to\{D' + 2D)^ 

\ mi mi J 


since |/q-(s)| < G^ (defined in Section III-A). We define A as 

A = AniaxG ai^Oiup Oilow)' 
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Now, for the third step of our proof we make 


Aoo{t) - ^oo(fo)| 


< 


< 


A{ 2- Tnz “h 2- \t — + 2Gj'\t — to\{D' + 2D) 

mi mi 


1 

m’ 


for any choice of m > 0. We do this hy making each of the three terms less than l/(3m). 
For the first term we want 

2AGTTnz ^ 
mi 3m 


We solve and get 


mi > QmAGxTnz- 


Since for any value of mi > 0 we can find an n > 0, this condition can he satisfied. 
For the third term we want 

2AGTiD'+ 2D)\t - to\ < 


This gives us 

~ ^ 6AGt{D^ + 2D)m' 

Since the only requirement is |f — fo| < 1/n for n chosen hy any given mi > 0, we can always choose 
\t — to\ as small as desired. Thus, this condition can he satisfied. 

With the second term we want the condition 


2AGt 

mi 




to\ 


1 

3m 


which means that 

\t - to\ ^ 1 

mi GmAGx 

Again, this condition can he satisfied since we can choose mi as large as we want and \t — to\ as small 
as we want as long as |f — to\ < 1/n for a given mi. 

Thus, for any m > 0, we first choose mi > QmAGxTnz- Then, we find an n' > 0 such that \t — to\ < 
Ijn' implies that \p{t) — p{to)\ < 1/mi for all t and to if we remove the discontinuities in p{t). Then, 
if necessary, n' is increased to n so that \t — to\ < l/n implies that \t — to\ < 1/{6AGt{D' + 2D)m) 
and \t — to\/mi < l/( 6 mAGr). If no increase is necessary, then n = n'. With this choice of n > 0, 
l^oo(f) — ^oo(fo)| < 1/m. As a result, for any m, we can find an n such that \t — to\ < Ijn implies 
that |Aoo(f) — Aoo(fo)| <'^lm. Thus, Aoo(f) is continuous. 
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This completes the proof for Lemma 3. A 

Proof of Theorem 2. Let us start hy writing (30) as 


N 


= E 


-^max Kfix ^j,i 

jr~ 


pit To Ti Dfix ^j,i) — 


1 


N 


i=l 2=1 

where Mi{t, s) = AmaxKfixKj^ip{t — To — Ti — Dfix — Djf). Recall that the dependence on s comes 
from the fact that the density of Tj is a function of a* which is characterized hy fa{s). This notation is 
analogous to the notation used in Section III-C. Following the steps in the proof of Lemma 1, we can 
quickly show that the limiting aggregate waveform at node j will take on the form 


pO'.up 

Vit)= E{Mi{t,s))fa{s)ds, 

J Otlo-w 

where 


(43) 


E{NU{t,s)) 

/ oo rO roo 

/ / g{-y)g{.x)p{t - to - V' - y - x)fD, s)dxdyd'tp, 

-oo J —oo J 0 

with g{-) = Therefore, we can prove Theorem 2 in two steps: 

• To show that r]{t) is odd about tq, we need to show that E{Mi(t, s)) is odd in t about tq, i.e. 

E{Mi{TQ + ^, s)) = -E{Mi{To -^,s)) for ^ > 0. 

• To show a zero-crossing at tq, show that E{Mi{To, s)) = 0. 

These two steps come directly from the form of r]{t) in (43). 

We first show that E{Mi{To + C,s)) = —E{Mi{TQ — ^,s)) for ^ > 0. Using the fact that Kfix = 
K{S~^{-Dfix)) = g{-Dfix) and Kjf = g{Dj^i), we have the following: 


E{Mi{To + ^,s)) 

— E(^Amaxgi~Dfix)g{Djf)p[^ — [Ti + Dfix T 

= ~E[Axnaxg{—Dfix)g{Dj^i)p{—^ + [Ti + Dfix + DjA)) 
OO /•O /•oo 


(A) ^ 


(c) 


/ OO /•u poo 

/ / g{-y)g{x)p{-i + [A+ y+ x\)fDAx)fDf,Ay)fT{A,s)dxdyd'ii) 

-oo J —oo J 0 

f 

J oo 

/ oo pO poo 

/ / g{-u)g{z)p{-i -[w + u + 2])/d, {u)fT{w, s)dzdudw 

-OO J —OO J 0 


—OO /•O /•— oo 

max! II g{z)g{-u)p{-^-[w + z + u])fDA-u)fDfiA-z)fTi-w,s)dudzdw 

oo J OO Jo 

OO /“O roo 


= -E[Araaxg{-Dfix)g{Dj^i)p{-i - [Ti + Dfix + DfA)) 
= -E{MAto-Ls)), 
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where (a) follows because p{t) = —p{—t) and at (6) we did a change of variables with u = —x, 
w = —'ijj, and z = —y. (c) follows from /t(x,s) = /t(— x,s) and fDj{x) = We thus have 

E{Mi{To + s)) = -E{Mi{TQ - s)) for ^ > 0 . 

E{Mi{TQ, s)) = 0 can now be shown as follows. Using the just proven fact that E{Mi{TQ + s)) = 

—E{Mi{To — ^,s)) for ^ > 0, setting ^ = 0 gives us E{Mi{To,s)) = —£'(Mj(ro, s)). This implies that 
E{Mi{To,s)) = 0 . 

This completes the proof for Theorem 2. A 
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